www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/fresnel.m

    function varargout = fresnel( U1,z,varargin)
%FRESNEL Fresnel diffraction
%  Syntax:
%  [U2,dxy2] = fresnel(U1,z,dxy1,lambda)
%  U2 = fresnel(U1,z)
%  fresnel(U1,z,...)
%
%  U1 is the wavefront of the object plane
%  U2 is the wavefront of the diffraction plane
%  U1 and U2 are all two-dimensional array
%  size of U1 and U2 are even
%  dxy1 is the sampling distance on the object plane. dxy1==[dx1,dy1]
%  if dxy1 is not input, it is set to [0.01,0.01]
%  dxy2 is the sampling distance on the diffraction plane. dxy2==[dxy2(1),dy2]
%  z is the distance between object plane and diffraction plane
%  z can be a scalar or a 2-element vector [zx,zy] to eliminate astigmatism
%  lambda is the wavelength of the laser
%
%  if there is no output, image on diffraction plane will be displayed
%  else, no image will be displayed
%
%  the origin of coordinates is at M/2+1,N/2+1
%
%  the relation between sampling distance on the two planes is:
%  dx2=lambda*z/dx1/N, dy2=lambda*z/dy1/M
%
error(nargchk(2,4,nargin))
if nargout>2
    error('Too many output arguments')
end

if length(z)==1
    z=[z,z];
end

[M,N]=size(U1);
switch nargin
    case 2
        dxy1=[0.01,0.01];
        lambda=5.32e-4;
    case 3
        dxy1=varargin{1};
        lambda=5.32e-4;
    case 4
        dxy1=varargin{1};
        lambda=varargin{2};
end

if isempty(dxy1)
    dxy1=[0.01,0.01];
end

k=2*pi/lambda;
dxy2=abs(lambda.*z./dxy1./[N,M]);

%-----------------------------------------------------
phase1=zeros(M,N);
for m=1:M
    for n=1:N
        x1=(n-N/2-1)*dxy1(1);
        y1=-(m-M/2-1)*dxy1(2);
        phase1(m,n)=exp(i*k/2*(x1^2./z(1)+y1^2./z(2)));
    end
end
U1=U1.*phase1;
clear phase1
%-----------------------------------------------------
phase2=zeros(M,N);
for m=1:M
    for n=1:N
        x2=(n-N/2-1)*dxy2(1);
        y2=-(m-M/2-1)*dxy2(2);
        phase2(m,n)=exp(i*k/2*(x2^2./z(1)+y2^2./z(2)));
    end
end
%-----------------------------------------------------
temp=zeros(M,N);
for m=1:M
    for n=1:N
        temp(m,n)=(-1).^(m+n);
    end
end
z=(z(1)+z(2))/2;
if z>0
    U2=exp(i*k*z)/(i*lambda*z)*dxy1(1)*dxy1(2)*(-1)^((M+N)/2)*phase2.*temp.*fft2(U1.*temp);
else
    U2=exp(i*k*z)/(i*lambda*z)*dxy1(1)*dxy1(2)*(-1)^((M+N)/2)*phase2.*temp.*ifft2(U1.*temp)*M*N;
end
switch nargout
    case 0
        u2=abs(U2);
        OpticImage(u2,dxy2(1),dxy2(2));xlabel('x');ylabel('y');title('Amplitude on Diffraction Plane');
    case 1
        varargout{1}=U2;
    case 2
        varargout{1}=U2;
        varargout{2}=dxy2;
end